
################################
### Regional income measures ###
################################

rm(list = ls())
gc()

library(data.table)
library(tidyverse)

### Set WD
setwd("~/Dropbox (Princeton)/Research/*Spatial_Inequality/paper_SMD_inequality/JoP_final/replication")  

### Regional income
load("data/inc_comb.RData")

inc_comb_sd <- inc_comb[ year >= 2000, lapply(.SD, function(x) sd(x, na.rm=T) ) , by=c("countryname","year"), .SDcols=c("income") ]
setnames(inc_comb_sd, "income","sd_income")

pdf("out/si_fig_d1.pdf", width=8, height=6.5) # sd_income_v1.pdf
(plot <- (ggplot(data=inc_comb_sd[ !is.na(sd_income) ], aes(y=log(sd_income), x=year))
          + geom_line()
          + theme_bw()
          + facet_wrap(~ countryname)
          + ylab("Log standard deviation of regional income")
          + xlab("Year")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), 
                  panel.grid = element_blank(), axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
          
))
dev.off()


### Annual change
inc_comb_sd_chg <- inc_comb_sd[ !is.na(sd_income) & !is.na(countryname)]

setkey(inc_comb_sd_chg, countryname, year)

inc_comb_sd_chg$countryname2 <- factor(inc_comb_sd_chg$countryname, levels=rev(unique(inc_comb_sd_chg$countryname)))

inc_comb_sd_chg[ , year := as.character(year) ]

# Calculate percentage change from one year to another
inc_comb_sd_chg[ , sd_income_lag_1yr := shift(.SD, type='lag', n=1), by = "countryname", .SDcols="sd_income"]
inc_comb_sd_chg[ , sd_income_pct_chg := (sd_income - sd_income_lag_1yr) / sd_income_lag_1yr ]

inc_comb_sd_chg$countryname2 <- factor(inc_comb_sd_chg$countryname, levels=rev(unique(inc_comb_sd_chg$countryname)))

pdf("out/si_fig_d2a.pdf", width=5, height=5)
(plot <- (ggplot(data=inc_comb_sd_chg, aes(y=sd_income_pct_chg, x=countryname2))
          + geom_hline(yintercept=0, colour = "gray50", linetype="dashed", size=.5)
          + geom_boxplot(fill="blue")
          + theme_bw()
          + coord_flip()
          + ylab("Annual change in regional inequality")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), 
                  panel.grid = element_blank(), axis.title.y = element_blank())
          
))
dev.off()


pdf("out/si_fig_d2b.pdf", width=5, height=5)
(plot <- (ggplot(data=inc_comb_sd_chg[ !is.na(sd_income) ], aes(y=log(sd_income), x=countryname2))
          + geom_boxplot(fill="blue")
          + theme_bw()
          + coord_flip()
          + ylab("Log regional income")
          + xlab("Year")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), 
                  panel.grid = element_blank(), axis.title.y = element_blank())
          
))
dev.off()



